module rhofunc

    implicit none
    public

    integer, parameter :: dp = kind(0.d0)

contains

    pure real(dp) function rho(x, y)
        real(dp), intent(in) :: x, y

        if (x > 0.6_dp .and. x < 0.8_dp .and. y > 0.6_dp .and. y < 0.8_dp) then
            rho = 1.0_dp
        else if (x > 0.2_dp .and. x < 0.4_dp .and. y > 0.2_dp .and. y < 0.4_dp) then
            rho = -1.0_dp
        else
            rho = 0.0_dp
        end if

    end function

end module rhofunc
